Recurrence Relations and Fast Algorithms 



We construct fast algorithms for evaluating transforms associ- 
ated with families of functions which satisfy recurrence relations. 
These include algorithms both for computing the coefficients in lin- 
ear combinations of the functions, given the values of these linear 
combinations at certain points, and, vice versa, for evaluating such 
linear combinations at those points, given the coefficients in the lin- 
ear combinations; such procedures are also known as analysis and 
synthesis of series of certain special functions. The algorithms of 
the present paper are efficient in the sense that their computational 
costs are proportional to n(lnn) (ln(l/e))^, where n is the amount 
of input and output data, and e is the precision of computations. 
Stated somewhat more precisely, we find a positive real number C 
such that, for any positive integer n > 10 and positive real num- 
ber e < 1/10, the algorithms require at most Cn (Inn) (ln(l/£))^ 
floating-point operations and words of memory to evaluate at n ap- 
propriately chosen points any linear combination of n special func- 
tions, given the coefficients in the linear combination, where e is the 
precision of computations. 

1 Introduction 

Over the past several decades, the Fast Fourier Transform (FFT) and its 
variants (see, for example, [11]) have had an enormous impact across the 
sciences. The FFT is an efficient algorithm for computing, for any positive 
integer n and complex numbers /3i, /32, • ■ ■ , Pn-i, Pn, the complex numbers 
ai, Q!2, . . . , a-a-i, ctn defined by 
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for J = 1, 2, . . . , n — 1, n, where /i, /2, . . . , fn-i, fn are the functions 
defined on [—1, 1] by 

fkix) = exp —j (2) 

for fc = 1, 2, . . . , n — 1, n, and Xi, X2, ■ ■ ■ , Xn-i, Xn are the real numbers 
defined by 



for fc = l, 2, — 1, n. The FFT is efficient in the sense that there 

exists a reasonably small positive real number C such that, for any positive 
integer n > 10, the FFT requires at most C n Inn floating-point operations 
and words of memory to compute ai, a2, • • • , ctn-i, ctn in (1) from (3i, (32, 
. . . , Pn-ii 0n- In contrast, evaluating the sum in (1) separately for every 
j = 1, 2, . . . , n — 1, n costs at least operations in total. 

The present paper introduces similarly efficient algorithms for comput- 
ing ai, a2, an in (1) from (32, (3n-i, (3n, and (when 
appropriate) for the inverse procedure of computing (3i, (32, ■ ■ ■ , Pn-i, /?« 
from ai, a2, • • • , an-i, a„, for more general collections of functions /i, /2, 
fn-i, fn and points xi, X2, Xn-i, Xn than those defined in (2) 
and (3). Specifically, the present paper constructs algorithms for classes 
of functions satisfying recurrence relations. The present paper describes in 
detail a few representative examples of such classes of functions, namely 
weighted orthonormal polynomials and Bessel functions of varying orders. 
These collections of functions satisfy recurrence relations of the form 

g{x) fk (x) = Ck-i fk-i{x) + dk fk (x) + Ck Ik+i{x) (4) 

for all X in the domain, where Ck-i, Ck, and dk arc real numbers and cither 
g{x) = X OT g{x) = i; Cfe, dk, and g vary with the collection of fimctions 
under consideration. 

The algorithms of the present paper all rely on the following two obser- 
vations: 

1. The solutions to the recurrence relation (4) are the eigenvectors cor- 
responding to eigenvalues g{x) of certain tridiagonal real self-adjoint 
matrices. 

2. There exist fast algorithms for determining and applying matrices 
whose columns are normalized eigenvectors of a tridiagonal real self- 
adjoint matrix, and for applying the adjoints of these matrices of 
eigenvectors. 



2 



The first observation has been well known to numerical analysts at least 
since the seminal [3] appeared; the second observation has been reasonably 
well known to numerical analysts since the appearance of the celebrated [5]. 
However, the combination seems to be new. 

The methods described in the present paper should lead to fairly effi- 
cient codes for computing a variety of what are known as (pseudo) spectral 
transforms. In particular, we can use the methods to construct fast al- 
gorithms for calculations involving spherical harmonics (see Remark 37 
below). 

We refer the reader to [13] and its compilation of references for prior 
work on related fast algorithms, as well as to [7] for an alternative approach 
that is suitable for certain applications, and to [9] for its refined accounting 
of computational costs. The present paper introduces techniques that are 
substantially more efficient than the extremely similar ones for which [13] 
reports on far-from-optimal implementations. We intend to report sepa- 
rately on carefully optimized implementations of the techniques described 
in the present paper, based in part upon the approach introduced in [8]. 
We gave a preliminary version of the present paper in [15]. 

The present paper has the following structure: Subsection 2.1 summa- 
rizes properties of fast algorithms for spectral representations of tridiagonal 
real self-adjoint matrices. Subsection 2.2 reiterates facts having to do with 
recurrence relations for orthonormal polynomials, Subsection 2.3 reiterates 
facts having to do with recurrence relations for Bessel functions, and Sec- 
tion 3 employs the subsections of Section 2 to construct fast algorithms for 
various purposes. 

2 Preliminaries 

This section summarizes certain widely known facts from numerical and 
mathematical analysis, used in Section 3. 

2.1 Divide-and-conquer spectral methods 

This subsection summarizes properties of fast algorithms introduced in [4] 
and [5] for spectral representations of tridiagonal real self-adjoint matrices. 
Specifically, there exists an algorithm such that, for any tridiagonal real 
self-adjoint matrix T, (firstly) the algorithm computes the eigenvalues of 
T, (secondly) the algorithm computes any eigenvector of T, (thirdly) the 
algorithm applies a square matrix U consisting of normalized eigenvectors 
of T to any arbitrary column vector, and (fourthly) the algorithm applies 
U"^ to any arbitrary cohimn vector, all using a number of operations and 
words of memory proportional to n (Inn) (ln(l/e))'', where n is the positive 
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integer for which T and U are nxn, and e is the precision of computations. 
The following is a more precise formulation. 

For any positive integer n, self-adjoint nxn matrix T, and real n x 1 
column vector v, we define ||Tl| to be the largest of the absolute values of 
the eigenvalues of T, 6t to be the minimum value of the distance |A — /Lt| 
between any two distinct eigenvalues A and of T, and 

n 

(5) 

k=l 

where vi, V2, ■ ■ ■ , Un-i, siie the entries of v; we say that v is normalized 
to mean that ||u|| = 1. As originated in [5], there exist an algorithm and a 
positive real number C such that, for any positive real number e < 1/10, 
positive integer n > 10, tridiagonal real self-adjoint nxn matrix T with n 
distinct eigenvalues, real unitary matrix U whose columns are n normalized 
eigenvectors of T, and real n x 1 column vector v, 

1. the algorithm computes to absolute precision ||T|| e the n eigenvalues 
of T, using at most 

Cn(lnn)(ln(l/£))^ (6) 
floating-point operations and words of memory, 

2. the algorithm computes to absolute precision ||r|| e/St the n en- 
tries of the matrix- vector product Uv, using at most 

Cn(lnn)(ln(l/e))=^ (7) 

operations and words of memory, 

3. the algorithm computes to absolute precision ||T|| ||u|| e/^r the n en- 
tries of the matrix-vector product U'^ v, using at most 

Cn(lnn)(ln(l/£))=^ (8) 

operations and words of memory, and, 

4. after the algorithm performs some precomputations which are par- 
ticular to T at a cost of at most 

Cn{\nn){ln{l/s)f (9) 

operations and words of memory, the algorithm computes to absolute 
precision ||T|| £/(5r the kn entries of any k normalized eigenvectors 
of T, using at most 

Ckn{ln{l/s)f (10) 
operations and words of memory, for any positive integer k. 
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Remark 1 Wc omitted distracting factors of very small powers of n in the 
precisions mentioned in the present subsection. Also, the bounds on the 
number of operations and words of memory are extremely conservative; in 
actual implementations the running-times of the algorithm appear to scale 
much better with respect to the precision e. 

Remark 2 In the second item of the present subsection, the algorithm in 
fact requires at most 

Cfcn(lnn)(ln(l/e))2 (11) 

operations and words of memory to compute the matrix-vector products 

U , Uv'^, . . . , U f*^"^, Uv'', for any positive integer fc, and real n x 1 
column vectors v^, , v^~^ , , after the algorithm performs some 
precomputations which are particular to T at a cost of at most 

Cn (Inn) (ln(l/e))-'^ (12) 

operations and words of memory. Moreover, we can improve the precisions 
to which the algorithm calculates U v^, Uv"^, . . . , Uv''~^, Uv'','by perform- 
ing more expensive precomputations (using higher-precision floating-point 
arithmetic or precomputation algorithms whose costs are not proportional 
to n Inn, for example). Similar considerations apply to the third item of 
the present subsection. 

Remark 3 There exist similar algorithms when the eigenvalues of T are 
not all distinct. 

2.2 Orthonormal polynomials 

This subsection discusses several classical facts concerning orthonormal 
polynomials. All of these facts follow trivially from results contained, for 

example, in [14]. 

Lemmas 8, 9, and 10, which formulate certain simple consequences of 
Theorems 4 and 7, are the principal tools used in Subsections 3.1 and 3.3. 
Lemmas 6 and 17 provide the results of some calculations for what are 
known as normalized Jacobi polynomials, a classical example of a family 
of orthonormal polynomials; the results of analogous calculations for some 
other classical families of polynomials are similar and therefore have been 
omitted. The remaining lemmas in the present subsection. Lemmas 12 
and 15, deal with certain conditioning issues surrounding the algorithms 
in Subsections 3.1 and 3.3 (see Remark 16). The remaining theorem in the 
present subsection. Theorem 14, describes what are known as Gauss-Jacobi 
quadrature formulae. 
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In the present subsection, we index vectors and matrices starting at 
entry 0. 

We say that a is an extended real number to mean that a is a real 
number, a = +00, or a = —00. For any real number a, we define the 
intervals [a, 00] = [a, 00) and [—00, a] = (—00, a]; we define [—00,00] = 
(—00, 00). 

For any extended real numbers a and b with a < b and nonncgative 
integer n, we say that po, pi, . . . , Pn-i, Pn (f^e orthonormal polynomials on 
[a, b] for a weight w to mean that «; is a real-valued nonnegative integrable 
function on [a, 6], pk is a polynomial of degree k, the coefficients of x", , 
. . . , x''~^, in Pk{x) are real, and the coefficient of x'' in Pk{x) is positive 
for fc = 0, 1, . . . , n — 1, n, and 

^ dx w{x) pj{x) pk{x) = I •j! ~ ^ (13) 

for j, fc = 0, 1, . . . , n — 1, n. 

The following theorem states that a system of orthonormal polynomials 
satisfies a certain three-term recurrence relation. 

Theorem 4 Suppose that a and b are extended real numbers with a < b, n 
is a positive integer, and po, pi, . . . , Pn-i, Pn are orthonormal polynomials 
on [a, b]. 

Then, there exist real numbers co, ci, . . . , c„_2, c„_i and do, di, . . . , 
dn-2, dn-i such that 

xpo{x) = doPo{x) + coPi{x) (14) 

for any x G [a,b], and 

xpk{x) = Ck-iPk-i{x) +dkPk{x) +CkPk+i{x) (15) 
for any x G [a, b] and fc = 1, 2, . . . , n — 2, n — 1. 

Proof. Theorem 3.2.1 in [14] provides an equivalent formulation of the 
present theorem. □ 

Remark 5 In fact, Cfc > for A; = 0, 1, . . . , n — 2, n — 1, in (14) and (15). 

The following lemma provides expressions for cq, Ci, c„_2, c„_i 
and do, rfi, dn-2, dn-i from (14) and (15) for what are known as 
normalized Jacobi polynomials. 



6 




Lemma 6 Suppose that a = — 1, 6=1, a and /? are real numbers with 
a > —1 and /3 > —1, n is a positive integer, and po, Pi, ■ . ■ , Pn-i, Pn are 
the orthonormal polynomials on [a, b] for the weight w defined by 

w{x) = {1 - x)" {1 + x)/^ . (16) 

Then, 

4(fc + l)(fc + a + l)(fc + /3+l)(fc + a + /3 + l) 
(2fc + a + /3 + l)(2fc + a + /3 + 2)2(2A; + a + /3 + 3) ^ ' 

* " (2A: + a + /3)(2fc + a + /3 + 2) ^^^^ 

for A; = 0, 1, . . . , n — 2, n — 1, where cq, Ci, . . . , c„_2, c„_i and do, di, . . . , 
d„_2, are from (14) and (15). 

Proof. Formulae 4.5.1 and 4.3.4 in [14] together provide an equivalent 
formulation of the present lemma. □ 

The following theorem states that the polynomial of degree n in a sys- 
tem of orthonormal polynomials on [a, b] has n distinct zeros in [a, b]. 

Theorem 7 Suppose that a and b are extended real numbers with a < b, n 
is a positive integer, and Pq, pi, . . . , Pn-i, Pn are orthonormal polynomials 
on [a, b]. 

Then, there exist distinct real numbers xq, Xi, Xn-2, Xn-i such 
that Xk G [a, b] and 

PniXk) = (19) 
for = 0, 1, . . . , n — 2, n — 1, and 

Xj Xk (20) 

when j ^ k iov j,k = 0, 1, . . . , n — 2, n — 1. 

Proof. Theorem 3.3.1 in [14] provides a slightly more general formulation 
of the present theorem. □ 

Suppose that a and b are extended real numbers with a < 6, n is a 
positive integer, and Po, Pi, . . . , Pn-i, Pn are orthonormal polynomials on 
[a,b] for a weight w. We define T to be the tridiagonal real self-adjoint 
n X n matrix with the entry 



Cj-1, k=j-l 

h — n 

(21) 



dj, k = j 



Cj, k = j + l 
0, otherwise (when fc < j — 1 or fc > j -|- 1) 
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for j,k = 0, 1, . . . , n — 2, n — 1, where cq, ci, . . . , c„_2, c„_i and do, di, 
. . . , dn-2, dn-1 are from (14) and (15). For k = 0, 1, . . . , n — 1, n, we 
define the function qi- on [a, b] by 

qk{x) = y/w(x) pk{x). (22) 
We define U to be the real n x n matrix with the entry 

Uj,k = -^ik^= (23) 



for j, fc = 0, 1, . . . , n — 2, n — 1, where go, ■ • ■ , <ln-2, Qn-i are defined 
in (22), and xq, xi, . . . , Xn-2, Xn-i are from (19). We define A to be the 
diagonal real n x n matrix with the entry 

for j,k = 0, 1, . . . , n — 2, n — 1, where xq, xi, . . . , Xn-2, Xn-i are from (19). 
We define S to be the diagonal real nx n matrix with the entry 



S.k = { \l^"rn=o{lm{Xj)f , k^J (25) 

I 0, k^j 

for j, fc = 0, 1, . . . , n — 2, n — 1, where qo, 9i, ■ • ■ , qn-2, Qn-i are defined 
in (22), and xq, xi, . . . , a;„_2: in-i are from (19). We define e to be the 
real nxl column vector with the entry 

^'==1 0, k^Q (26) 

for fc = 0, 1, . . . , n — 2, n — 1. 

The following lemma states that J7 is a matrix of normalized eigenvec- 
tors of the tridiagonal real self-adjoint matrix T , and that A is a diagonal 
matrix whose diagonal entries are the eigenvalues of T (which, according 
to (20), are distinct). 

Lemma 8 Suppose that a and b are extended real numbers with a < b, n 
is a positive integer, and Po, Pi, ■ ■ , Pn-i, Pn are orthonormal polynomials 
on [a, b] for a weight w. 
Then, 

U'^TU = K, (27) 

where T is defined in (21), U is defined in (23), and A is defined in (24). 
Moreover, U is real and unitary. 
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Proof. Combining (14), (15), and (19) yields that 

TU = UK. (28) 

Combining (28), (23), (24), and (20) yields that f/ is a real matrix of nor- 
malized eigenvectors of T, with distinct corresponding eigenvalues. There- 
fore, since eigenvectors corresponding to distinct eigenvalues of a real self- 
adjoint matrix are orthogonal, U is orthogonal. Applying U'^ from the left 
to both sides of (28) yields (27). □ 

The following lemma expresses in matrix notation the analysis and 
synthesis of linear combinations of weighted orthonormal polynomials for 
which Subsection 3.3 describes fast algorithms. 

Lemma 9 Suppose that a and b are extended real numbers with a < b, n 
is a positive integer, po, Pi, • • • , Pn-i, Pn are orthonormal polynomials on 
[a, b] for a weight w, and a and (3 are real n x 1 column vectors, such that 
a has the entry 

n-l 

aj = Y.(^kqk{xj) (29) 

fc=0 

for J = 0, 1, . . . , n — 2, n — 1, where Qq, qi, . . . , qn-2, Qn-i are defined 
in (22), and Xo, Xi, . . . , Xn-2, Xn-i are from (19). 
Then, 

a = SU'^p (30) 

and 

/3 = C/5-^Q!, (31) 

where U is defined in (23), S is defined in (25), and and U S-'^ a 

are matrix-matrix- vector products. 

Proof. Combining (23) and (25) yields (30). According to Lemma 8, U 
is real and unitary. Therefore, applying the matrix-matrix product U 
from the left to both sides of (30) yields (31). □ 

The following two lemmas provide alternative expressions for the entries 
of 5 defined in (25). 

Lemma 10 Suppose that a and b are extended real numbers with a < b,n 
is a positive integer, and po, pi, . . . , Pn-i, Pn are orthonormal polynomials 
on [a, b] for a weight w. 
Then, 

Sk,k = ^ (32) 



(C/T e)k !l dx w{x) 
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for fc = 0, 1, . . . , n — 2, n — 1, where S is defined in (25), U is defined 
in (23), e is defined in (26), {U^ e)o, {U^ e)i, . . . , (t/'^ e)„_2, {W^ e)„_i are 
the entries of the matrix- vector product U'^ e, and xo, xi, . . . , Xn-2, Xn-i 
are from (19). 

Proof. Combining (23) and (26) yields that 

(U e)k = , _^ == (33) 
yTZ=o ilmiXk)) 

for A; = 0, 1, . . . , n — 2, n — 1. Since the polynomial po has degree 0, 
combining (13) and (22) yields that 



.o(.) = (34) 

V la dy My) 

for any x € [o, b]. Combining (25), (33), and (34) yields (32). □ 



Remark 11 Formula 2.6 in [3] motivated us to employ the equivalent (32). 

Lemma 12 Suppose that a and b are extended real numbers with a < b, 
n is a positive integer, po, Pi, • • • , Pn-i, Pn arc orthonormal polynomials 
on [a,b] for a weight w, and fc is a nonnegative integer, such that Inw is 
differentiable at the point Xk from (19). 
Then, 

{Sk,kf = Cn-l qn-l{Xk) '^1n{Xk), (35) 

where Sk,k is defined in (25), c„_i is from (15), g„_i and Qn are defined 
in (22), and Xk is from (19). 

Proof. Formula 3.2.4 in [14] provides a slightly more general formulation 
of the present lemma. □ 



Remark 13 There exist similar formulations of Lemma 12 when it is not 
the case that Inw is differentiable at Xk- 

The following theorem describes what are known as Gauss- Jacobi quad- 
rature formulae for orthonormal polynomials. 
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Theorem 14 Suppose that a and h are extended real numbers with a < b, 
n is a positive integer, and j30) Pit ■ ■ • , Pn-iiPn are orthonormal polynomials 
on [a, b] for a weight w. 

Then, there exist positive real numbers wq, wi, . . . , Wn-it called 

the Christoffel numbers for xq, Xi, . . . , Xn-2, Xn-i, such that 

/ dx w{x) p{x) = ^ Wkp{xk) (36) 
fe=o 

for any polynomial p of degree at most 2n— 1, where Xo,Xi, . . . , a;„_2, a;„_i 
arc from (19). 

Proof. Theorems 3.4.1 and 3.4.2 in [14] together provide a slightly more 
general formulation of the present theorem. □ 

The following lemma provides alternative expressions for the entries of 
S defined in (25). 

Lemma 15 Suppose that a and b are extended real numbers with a < b, n 
is a positive integer, and po, Pi, • • • , Pn-i, Pn are orthonormal polynomials 
on [a, b] for a weight w. 
Then, 

{Sk^kf = ^ (37) 

Wk 

for fc = 0, 1, . . . , n — 2, n — 1, where S is defined in (25), and wq, Wi, . . . , 
Wn-2, Wn-i are the Christoffel numbers from (36) for the corresponding 

points xq, xi, . . . , x„_2, Xn-i from (19). Moreover, there exist extended 
real numbers yo, yi, yn-i, Vn such that a = yo < yi < ■ ■ ■ < yn-i < 
yn = b and 

rVk+i 

'Vk 

for fc = 0, 1, . . . , n — 2, n — 1. 



rVk+i 

I dxw{x) (38) 

J Vk 



Proof. Formula 3.4.8 in [14] provides an equivalent formulation of (37). 
Formula 3.41.1 in [14] provides a slightly more general formulation of (38). 
□ 



Remark 16 The formulae (25), (35), (37), and (38) give some insight into 
the condition number of S. For instance, due to (25), the entries of S are 
usually not too large. 

The following lemma provides an alternative expression for the entries 
of S defined in (25) for what are known as normalized Jacobi polynomials. 
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Lemma 17 Suppose that a = —1, b = 1, a and /3 are real numbers with 
a > —1 and /? > —1, n is a positive integer, and po, pi, . . . , Pn-i, Pn are 
the orthonormal polynomials on [a, b] for the weight w defined by 



w{x) = {l-x)"{l+xf. 



(39) 



Then, 



Sk,k = 



1 



2n + a + /? + 1 



dx 



QniXk) 



(40) 



for A; = 0, 1, . . . , n — 2, n — 1, where S is defined in (25), xq, xi, . . . , 
from (19), and qn is defined in (22). 

Proof. Together with (37), Formulae 15.3.1 and 4.3.4 in [14] provide an 
equivalent formulation of the present lemma. □ 



2.3 Bessel functions 

This subsection disciisses several well-known facts concerning Bessel func- 
tions. All of these facts follow trivially from results contained, for example, 
in [17] and [10]. 

Lemmas 25, 26, 27, and 28 are the principal tools used in Subsec- 
tions 3.2 and 3.4. These lemmas formulate certain simple consequences of 
Theorems 19 and 22 and Corollary 23, by way of Lemmas 20 and 24. Lem- 
mas 32 and 33 provide closed-form results of some calciilations for what are 
known as spherical Bessel functions, a family of Bessel functions frequently 
encountered in applications. The remaining lemmas in the present subsec- 
tion, Lemmas 34 and 35, provide closed-form results of similar calculations 
for Bessel functions of arbitrary nonnegative orders (however, the results 
in Lemmas 34 and 35 arise from identities that are presumably not quite 
as familiar). 

In the present subsection, we index vectors and matrices starting at 
entry 1. 

Suppose that v is a. nonnegative real number. For any nonnegative 
integer k, we define the function fk on (0, oo) by 

fk{x) = ^ ^ Jt.+k{x), (41) 

x'^ 

where T is the gamma (factorial) function and J^+k is the Bessel function 
of the first kind of order u + k (see, for example, [17]). 
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Remark 18 Formula 8 of Section 3.1 in [17] provides a more general for- 
mulation of the fact that 

lim + J (3.) = 1 (42) 

which motivated our choice of normalization in (41). 

The following theorem states that /i, /2, /a, ... defined in (41) satisfy 
a certain three-term recurrence relation. 

Theorem 19 Suppose that v \s a. nonnegative real number. 
Then, 

- /i (x) = , 2''T{v + l) 1 ^ 

for any positive real number a;, and 

- fk{x) = ^ . , ^ fk-i{x) + ^ I, , /fe+i(a;) 

X 2i/(z/-Ffc- l)(i/-|-fc) 2^{y^k){y^k^\) 

(44) 

for any positive real number x and fc = 2, 3, 4, . . . , where /i, /2, /a, ... are 
defined in (41), F is the gamma (factorial) function, and is the Bessel 
function of the first kind of order v (see, for example, [17]). 

Proof. Formula 1 of Section 3.2 in [17] provides a somewhat more general 
formulation of the present theorem. □ 

Suppose that is a nonnegative real number and n is a positive integer. 
We define T to be the tridiagonal real self-adjoint n x n matrix with the 
entry 



2^(v+j-l){u+j) ' 

= , k=j + l 



, 0, otherwise (fc < j — 1, k — j, oy k > j + 1) 

(45) 

for J, fc = 1, 2, . . . , n — 1, n. For any positive real number x, we define 
V = v{x) to be the real n x 1 column vector with the entry 

Vk = , ^"^''^ (46) 



\/Em=l ifmix)f 
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for fc = 1, 2, . . . , n — 1, n, where /i, /2, . . . , /n-i, /n are defined in (41). 
For any positive real number x, we define 5 = 5{x) to be the real number 



2V(^' + n)(. + n+l) 

where /i, /2, . . . , /«, fn+i are defined in (41). 

The following lemma states that v is nearly an eigenvector of the tridiag- 
onal real self-adjoint matrix T corresponding to an approximate eigenvalue 
of I for any positive real number x such that Jv{x) = and S is small. 

Lemma 20 Suppose that u is a. nonnegative real number and n is a posi- 
tive integer. 



Tlien, 



{Tv)n Vr, 

X 



< S (48) 



and 



{Tv)k = -Vk (49) 



for k = 1, 2, . . . , n — 2, n — 1 and any positive real number x with 

J,{x) = 0, (50) 

where T is defined in (45), v = v{x) is defined in (46), {Tv)i, {Tv)2, ■ ■ ■ , 
{Tv)n-i, {Tv)n are the entries of the matrix-vector product Tv, S = 6{x) 
is defined in (47), and J^, is the Bessel function of the first kind of order v 

(see, for example, [17]). 

Proof. Combining (44), (43), and (50) yields (48) and (49). □ 



Remark 21 It is well known that, for any positive real number x, the 
quantity Jy+„+i(a;) and thence i5 defined in (47) decays extremely rapidly 
as n increases past a band around n = x of width proportional to x^^^; 
see, for example. Lemma 2.5 in [12], Chapters 9 and 10 in [1], or Chapter 8 
in [17]. Therefore, S is often small for x such that x <n and Jv{x) = 0. 

The following theorem states a simple Sturm sequence property of the 
eigenvalues of real self-adjoint tridiagonal matrices whose entries on the 
sub- and super-diagonals are nonzero. 

Theorem 22 Suppose that n is a positive integer, and T is a tridiagonal 
real self-adjoint n x n matrix, such that all entries on the sub- and super- 
diagonals of T are nonzero. 

Then, every eigenvalue of T has multiplicity 1. 



14 



Proof. Formula 7-7-1 in [10] provides an equivalent formulation of the 

present theorem. □ 



The following corollary of Theorem 22 states a simple Sturm sequence 
property of the eigenvalues of T defined in (45). 

Corollary 23 Suppose that is a nonnegative real number and n is a 
positive integer. 

Then, every eigenvalue of T defined in (45) has multiplicity 1. 

The following lemma bounds the distance between an approximate 
eigenvalue and the actual eigenvalue nearest to the approximation, as 
well as the discrepancy between the corresponding normalized approximate 
eigenvector and a corresponding normalized actual eigenvector. 

Lemma 24 Suppose that 7, A, and /x are real numbers, n is a positive 

integer, T is a real self-adjoint n x n matrix, and u and v arc real n X 1 
column vectors, such that A is the eigenvalue of T nearest to 7, A has 
multiplicity 1, /i is the eigenvalue of T nearest but not equal to A, 

Tu = Xu, (51) 

n 

Y^^Ukf = 1, (52) 



k=l 



and 



T^ivkf = 1. (53) 



fe=i 



Then, 



2 

|7-A|< E((^^)fe-7^fc) > (54) 

\ k=l 



and either (or both) 



for fc = 1, 2, . . . , n — 1, n, or 



Vk-Uk\< : -. (55) 

l/X — A| 



\-Vk-Uk\< . -. (56) 

l/i - A| 

for fc = 1, 2, . . . , n - 1, n, where {Tv)i, {Tv)^, {Tv)n-i, {Tv)n are 
the entries of the matrix- vector product Tv. 
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Proof. Formula 4-5-1 in [10] provides an equivalent formulation of (54). 

The proof of Formula 11-7-1 in [10], specifically Formula 11-7-3 in [10] 
and the comment immediately following Formula 11-7-3 in [10], provides a 
slightly more general formulation of the fact that (55) holds for k = 1, 2, 
. . . , n — 1, n, or that (56) holds for fc = 1, 2, . . . , n — 1, n. □ 

The following lemma bounds the changes in the eigenvalues and eigen- 
vectors induced by using the truncated matrix T defined in (45) (T is only 
n X n, not infinite-dimensional). 

Lemma 25 Suppose that A, fi, u, and x are real numbers, n is a positive 
integer, and w is a real n x 1 column vector, such that v > 0, x > 0, (50) 
holds, A is the eigenvalue of T nearest to i, /U is the eigenvalue of T nearest 
but not equal to A, 

Tu = Xu, 



and 



where T is defined in (45). 
Then, 



<s, 



and either (or both) 



\vk - Uk\ < 



25 



for fc = 1, 2, . . . , n — 1, n, or 



-Vk - Uk\ < 



25 



(57) 
(58) 

(59) 
(60) 

(61) 



for fc = 1, 2, . . . , n — 1, n, where v = v{x) is defined in (46) and 5 = 5{x) 
is defined in (47). 

Proof. Combining Corollary 23, (54), (48), and (49) yields (59). 

Combining Corollary 23, (55), (56), (48), and (49) yields that (60) holds 
for fc = 1, 2, . . . , n — 1, n, or that (61) holds for fc = 1, 2, . . . , n — 1, n. □ 

Suppose that v isa, nonnegative real number and n is a positive integer. 
We define a;i, 2:2, a;3, . . . to be all of the positive real numbers such that 



J^{xk) = 



(62) 
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for any positive integer k, ordered so that 

< < a;2 < < . . . , (63) 

where Ji, is the Bessel function of the first kind of order v (see, for ex- 
ample, [17]). We define S to be the diagonal real n x n matrix with the 
entry 




for j, fc = 1, 2, . . . , n - 1, n, where /i, /2, • • • , /„ are defined in (41), 

and xi, X2, ■ ■ ■ , a;„_i, x„ are defined in (62) and (63). We define e to be 
the real n x 1 column vector with the entry 

for fc = 1, 2, . . . , n — 1, n. 

The following lemma expresses in matrix notation the evaluations of 
linear combinations of Bessel functions for which Subsection 3.4 describes 
fast algorithms. 



Lemma 26 Suppose that v is & nonncgativc real number, n is a positive 
integer, and a and (3 are real n x 1 column vectors, such that a has the 
entry 

n 

'^i = ^l^kfk{xj) (66) 
fe=i 

for j = 1, 2, . . . , n — 1, n, where /i, /2, . . . , /n-i, fn are defined in (41), 
and xi, X2, ■ ■ ■ , Xn-i, x„ are defined in (62) and (63). 
Then, 

\a,-{SU^Ph\< ^f'^' (67) 
for any k — 1, 2, . . . , n — 1, n such that 

2 Sk,k S{xk) 



l/ifc — Afe| 



< \fi{xk)\, (68) 



where Xk is the eigenvalue of T defined in (45) nearest to — , Ufe is the 
eigenvalue of T nearest but not equal to Xk, S = S{xk) is defined in (47), 
?7 is a real n x n matrix whose fc'^ column is the normalized eigenvector 
of T corresponding to the eigenvalue Xk whose first entry has the same 
sign as fi{xk), S is defined in (64), and {SU'^(3)k is the k^^ entry of the 
matrix-matrix- vector product S U'^p. 
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Proof. Combining (60), (46), and (64) yields (67). 



□ 



The following two lemmas provide alternative expressions for the entries 
of 5 defined in (64). 



Lemma 27 Suppose that v is a, nonnegative real number and n is a posi- 
tive integer. 
Then, 



fl{Xk) 



< 



2 Sk,k 5{xk) 
\^J■k - Afel \{U^ e)k\ 



(69) 



for any fc = 1, 2, . . . , n — 1, n such that (68) holds, where S is defined 
in (64), /i is defined in (41), xi, X2, Xn-i, Xn are defined in (62) 
and (63), Afc is the eigenvalue of T defined in (45) nearest to fJ-k is the 
eigenvalue of T nearest but not equal to Xk, S = S{xk) is defined in (47), 
[/ is a real n x n matrix whose k^^ column is the normalized eigenvector 
of T corresponding to the eigenvalue whose first entry has the same 
sign as fi{xk), e is defined in (65), and {W^ e)k is the fc*^ entry of the 
matrix-vector product e. 



Proof. Combining (60), (46), and (65) yields that 



./i(.':A:: 



\/Em=l ifmiXk)f 



< 



2 Sjxk) 



(70) 



for any fc = 1, 2, . . . , n—1, n such that (68) holds. Combining (64) and (70) 
yields (69). □ 



Lemma 28 Suppose that u is a nonnegative real number and n is a posi- 
tive integer. 
Then, 

Mxk) = '-^^ -Mxk) (71) 

for fc = 1. 2, . . . . n — 1, n. where /i is defined in (41). xi, X2, • . • , .t„-i, .t„ 
are defined in (62) and (63), F is the gamma (factorial) function, and is 
the Bessel function of the first kind of order v (see, for example, [17]). 

Proof. Formula 4 of Section 3.2 in [17] provides a somewhat more general 
formulation of (71). □ 
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Remark 29 The right hand side of (69) involves the potentially trouble- 
some 

^ (72) 



However, due to (60), (46), and (64), if 

2 S{xk) 



(73) 



(74) 



l^fe — Afe| 

is small, then (72) is accordingly close to 

Sk,k 

which should not be unreasonably large. 

Remark 30 Numerical experiments indicate that — A^l in (67) and (69) 
is never exceedingly small for practical ranges of n; this is probably fairly 
easy to prove, perhaps using the properties of Sturm sequences. The fol- 
lowing remark appears to be relevant. 

Remark 31 Suppose that 1^ = and n is a positive integer. For any 
positive real number x, we define v = v{x) to be the real n x 1 column 
vector with the entry 

Vk = {-!)'' Vk (75) 

for fc = 1, 2, . . . , n — 1, n, where v = v{x) is defined in (46). Then, 
Formula 1 of Section 3.2 in [17] and Formula 2 of Section 2.1 in [17] lead 
to 

{Tv)r, + -Vn <S (76) 
X 

in place of (48), 

{Tv)k = --Vk (77) 

X 

for fc = 1, 2, . . . , n — 2, n — 1 in place of (49), etc. 

The following lemma states a special case of the Gegenbauer addition 
formula for Bessel functions. 

Lemma 32 Suppose that f = 
Then, 



E (/™(^))' = 2 (^^^ 

m=0 

for any positive real number x, where /o, /i, /2, ... are defined in (41). 
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Proof. Formula 3 of Section 11.4 in [17] provides a somewhat more general 
formulation of (78). □ 

The following lemma provides alternative expressions for the entries of 
S defined in (64) for what are known as spherical Bessel functions of the 
first kind. 

Lemma 33 Suppose that = i, e is a positive real number, and k and n 
are positive integers, such that 

oo 

E ifm{Xk)f<S, (79) 
m— n+1 

where fn+i, fn+2, fn+3, ■ ■ ■ are defined in (41), and Xk is defined in (62) 
and (63). 
Then, 



< e, (80) 



where Sk,k is defined in (64). 

Proof. Combining (64), (78), (79), and (62) yields (80). □ 

The following lemma provides a closed-form expression for the sum 
in (78), for any nonnegative order v. 

Lemma 34 Suppose that v is a nonnegative real number. 
Then, 

m=l ^ ' ^ ' 



(2z/+l)a; /2^r(i/ + l) 



Mx)hix) (81) 



for any positive real number x, where /i, /2, /s, ... are defined in (41), T 
is the gamma (factorial) function, and is the Bessel function of the first 
kind of order v (see, for example, [17]). 

Proof. Formula 57.21.1 of [6] provides a more general formulation of the 
present lemma. See also Formula 24 of [16] and the surrounding discussion 
for a self-contained derivation of the present lemma. □ 

The following lemma provides alternative expressions for the entries of 
S defined in (64). 
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Lemma 35 Suppose that v and e arc real numbers, and k and n are 
positive integers, such that i/ > 0, £ > 0, and (79) holds, where in (79), 
in+Xi fn+2, fn+3, ■ ■ ■ are defined in (41), and Xk is defined in (62) and (63). 
Then, 



{Sk,k) - 



2{v + l) 

where Sk,k is defined in (64), and /i is defined in (41). 
Proof. Combining (64), (81), (79), and (62) yields (82). 



(82) 



□ 



Remark 36 As in Remark 21, it is often possible to have e in (79), (80), 
and (82) be small for k such that Xk < n. 

3 Fast algorithms 

This section constructs efficient algorithms for computing the quadrature 
nodes and Christoffel numbers associated with orthonormal polynomials, 
for computing the zeros of Bessel functions, for the analysis and synthe- 
sis of linear combinations of weighted orthonormal polynomials, and for 
evaluations of linear combinations of Bessel functions. We describe the 
algorithms in Subsections 3.1 and 3.2 solely to illustrate the generality of 
the techniques discussed in the present paper; we would expect spcicial- 
ized schemes such as those in [2] to outperform the algorithms described in 
Subsections 3.1 and 3.2 in most, if not all, practical circumstances. Each 
subsection in the present section relies on both Subsection 2.1 and either 
Subsection 2.2 or Subsection 2.3. 

3.1 Quadrature nodes and Christoffel numbers associ- 
ated with orthonormal polynomials 

The entries of A in (27) are the nodes Xq, Xi, . . . , x„_2, a;„_i in (36). We 
can compute rapidly the entries of A in (27) using an algorithm as in the 
first item in Subsection 2.1, due to (27), since T in (27) is tridiagonal, real, 
and self-adjoint, U in (27) is real and unitary, and A in (27) is diagonal, 
with diagonal entries that according to (20) are distinct. For the same 
reason, we can apply rapidly the matrix U'^ to the vector e in (32) using 
an algorithm as in the third item in Subsection 2.1. We can then compute 
the Christoffel numbers wo,wi, . . . , Wn-2, Wn-i in (36) using (37) and (32). 
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3.2 Zeros of Bessel functions 



Wc can compute rapidly the zeros xi, X2, . . . , x^-i, a;„ defined in (62) 
and (63) for wfiicfi 6 defined in (47) is sufficiently small, using (59) and an 
algorithm as in the first item in Subsection 2.1, since T defined in (45) is 
tridiagonal, real, and self-adjoint, and (according to Corollary 23) has n 
distinct eigenvalues. 

3.3 Analysis and synthesis of linear combinations of 
weighted orthonormal polynomials 

Wc can apply rapidly the matrices U and U"^ in (31) and (30) using an 
algorithm as in the second and third items in Subsection 2.1, due to (27), 
since T in (27) is tridiagonal, real, and self-adjoint, U in (27) is real and 
unitary, and A in (27) is diagonal, with diagonal entries that according 
to (20) are distinct. Furthermore, we can apply rapidly the remaining 
matrices S and 5"^ in (30) and (31), since S and 5"^ in (30) and (31) 
are diagonal, once we use (37) and the algorithms from Subsection 3.1 to 
compute the entries of S and . 

Remark 37 Using the algorithms described in the present subsection, we 
can construct fast algorithms both for computing the coefficients in lin- 
ear combinations of spherical harmonics, given the values of these linear 
combinations at certain points, and, vice versa, for evaluating such linear 
combinations at those points, given the coefHcients in the linear combina- 
tions. We can handle spherical harmonics by constructing fast algorithms 
for what arc known as associated Legcndrc functions: sec, for example, [13]. 
For any nonnegative integers / and to, the normalized associated Legendre 
function of order m and degree I (often denoted by P; ) is equal to the 
function qi-m defined in (22) for the orthonormal polynomials on [—1,1] 
for the weight w defined by 

w{x) = {I - x)'^ {I + x)'^ . (83) 

Thus, we could utilize the algorithms discussed in the present subsection 
exactly as described, with (40) guaranteeing that the condition number 
of S is never too large for practical ranges of I and to. However, we 
would want to take; advantage of the symmetries of associated Legendre 
functions, by handling the even and odd functions separately, using the 
recurrence relation associated with P; [x) instead of the recurrence re- 
lation (15), which is associated with x Pi (./■). Wc; might also want to use 
the Christoffel-Darboux identity to compute interpolations to and from 
values at the zeros of various polynomials, as originated in [7] and [18], and 
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subsequently optimized and extended (see, for example, [8] and Remarks 11 
and 15 in [16]). 

3.4 Evaluations of linear combinations of Bessel func- 
tions 

We can apply rapidly the matrix U'^ in (67) using an algorithm as in the 
third item in Subsection 2.1, since T defined in (45) is tridiagonal, real, and 
self-adjoint, and (according to Corollary 23) has n distinct eigenvalues, and 
hence U in (67) can be chosen to be real and unitary. Furthermore, wc can 
apply rapidly the remaining matrix S in (67), since S in (67) is diagonal, 
once we use (69), an algorithm as in the third section in Subsection 2.1, 
and the algorithm from Subsection 3.2 to compute the entries of S. 
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